%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CS221 Programming Assignment 2
%%   Olga Russakovsky, Oct. 2009
%%   Stanford University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Uses entropy to decide which pixel to split on
%
% Input:
% - DigitSet is the set of all pixels at this node in the tree
%       DigitSet.pixels is a M x N matrix of pixel values, where M is the
%       number of example images and N is the number of pixels within each
%       example image 
%       DigitSet.labels is a M-dimensional vector of example labels (0-9)
%       DigitSet.weights is a M-dimensional vector of example weights
% - positiveLabel is the label of the positive class for this decision tree
% (0-9)
%
% Output:
% - pixelNum is the pixel number to split on (between 1 and N)
% - threshold is the split threshold (0.1, 0.2, ..., 0.9)


function [featureNum threshold] = choosePixelAndThreshold(DataSet)
   M = length(DataSet.y);  % number of examples
   if (M==0)
       error('choosePixelAndThreshold : data set is empty');
   end
   N = length(DataSet.x(1,:)); % number of features
  allWeights = sum(DataSet.weights);
  
  thresholds = DataSet.thresholds; % size T*N : 1 list of thresholds per feature
  T = DataSet.numThresholds;
       
  
       % create matrix of comparisons between pix>thresh
       a = reshape(DataSet.x,[M 1 N]);
       A = repmat(a, [1 T 1]); 
      % size(DigitSet.pixels)
       b = reshape(thresholds,[1 T N]);
       B = repmat(b, [M 1 1]);
       greaterThans = (A > B);
       % create total weights arrays
       totalweighted_bigger_than_thresh = repmat(DataSet.weights, [1 T N]) .* greaterThans; % dim M*T*N
       positiveweights = (DataSet.y == ones(M, 1)) .* DataSet.weights;  % dim M*1
       positiveweighted_bigger_than_thresh = repmat (positiveweights, [1 T N]) .* greaterThans;  % dim M*T*N
       m1 = sum(totalweighted_bigger_than_thresh, 1); % dim 1*T*N
       m1 = reshape(m1, [T N]);
       mplus1 = sum(positiveweighted_bigger_than_thresh, 1); % dim 1*T*N
       mplus1 = reshape(mplus1,[T N]);
       m2 = allWeights*ones(T,N) - m1;  % dim T*N
       mplus2 = sum(positiveweights)*ones(T,N) - mplus1; % dim T*N
       
       % compute probabilities, handling the case wgere m = 0
       zerom1 = (m1==zeros(T,N));
       zerom2 = (m2==zeros(T,N));
       p1 = mplus1 ./ (m1+zerom1) - mplus1 .* zerom1; % dim T*N : p=0 if m=0
       p2 = mplus2 ./ (m2+zerom2) - mplus2 .* zerom2; % dim T*N
       infogain = -m1 .* H_2D(p1) - m2 .* H_2D(p2); % dim T*N
       [maxVect, indexVect] = max(infogain); 
       [maxVal, indexVal] = max(maxVect);
       maxfeat = -1;
       maxthresh = -1;
       maxthresh = DataSet.thresholds(indexVect(indexVal),indexVal);
       maxfeat = indexVal;
       
       featureNum = maxfeat;
       threshold = maxthresh;
       
       if(maxfeat==-1)
           error('choosePixelAndThreshold : no splitting point found');
           maxfeat
           maxthresh
       else
           %disp('choosePixelAndThreshold : normal output, splitting point found');
       end
% CS221 TODO: implement this function
% -----------------------------------
% Debugging:
%   The following lines of code
%       DigitSet = loadDigits('data/training-10k-images.idx3', 'data/training-10k-labels.idx1');
%       [pixelNum threshold] = choosePixelAndThreshold(DigitSet, 3)
%   should return
%       pixelNum = 34
%       threshold = 0.2000


end


% Important note:
%  This function will get called millions of times as you are running
%  experiments; thus, you might want to consider some optimization tricks.
%  Be careful, start by first implementing the simplest possible solution
%  you can think of, and verify that your results are correct. 
%
% After you're done, you might want to consider that for-loops in Matlab
% are extremely inefficient, and so eliminating even just one inner loop
% will often greatly speed up your code. For instance, consider this trick: 
%    For each pixel p, you will be iterating over the possible thresholds
%    from 0.1 to 0.9 and computing the entory at each one. Letting
%           A = repmat(DigitSet.pixels(:, p), 1, length(thresholds))
%           B = repmat(thresholds, length(DigitSet.pixels(:, p)), 1)
%           greaterThans = (A > B);
%     will create a NxT array, where N is the number of training examples
%     and T is the number of thresholds, and greaterThans(i, j) will be 1
%     if the p^th pixel of training example i is greater than the j^th
%     threshold, and 0 otherwise.
%     You can then use it to construct other useful NxT arrays, e.g.,
%           weighted = repmat(DigitSet.weights, 1, length(thresholds)) .* greaterThans;
%     where weighted(i, j) will correspond to the weight of the i^th
%     training example if its p^th pixel value is greater than
%     thresholds(j), and 0 otherwise. You can then sum over the appropriate
%     dimension using
%           sum(weighted, 1)
%     to obtain an array of weights, one per threshold, that is useful for
%     this problem.
% You will need more than the arrays discussed here to efficiently
% implement the function, but this should help you get started. In general,
% you should:
%      (1) avoid using for-loops whenever possible and instead use
%          operations on arrays (e.g., 
%              [m i] = min(v)
%          will return the mininum element m of a vector v, and the
%          corresponding index i)
%      (2) avoid using if-statements if possible
%
% Please come by office hours if this is unclear or you need more help.

